#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include <time.h>
#include <sys/time.h>

#ifdef DEBUG
#define ERROR(fmt, args ...)	\
	printf("[Error] %s:%d\t" fmt, __FUNCTION__, __LINE__, ##args);
#endif	// DEBUG

#define TILE_WIDTH	16

static int SIZE;			// Matrix Size
static int ROW_SIZE;
float dL2Norm;
float* gMatrixA;
float* gMatrixB;
float* gMatrixResult;		// For test
float* gTargetRow;			// For test
float* gCoeffi; 			// For test

float* dTargetRow;			// For test
float* dCoeffi; 			// For test
float* dMatrixA;
float* dMatrixResult;		// For test



struct option input_option[] = { 
	{"size", 1, 0, 's'},
	{0, 0, 0, 0}
};
struct timeval startTime, endTime;

void show_matrix(float* gMatrix)
{
	int i, j;
	for(i=0; i<SIZE; ++i){
		for(j=0; j<SIZE+1; ++j){
			printf("%f\t", gMatrix[i*(SIZE+1)+j]);
		}
		printf("\n");
	}
	printf("\n");
}

int init()
{
	int i, j;

	// Make Input matrix
	gMatrixA = (float*)malloc(sizeof(float)*SIZE*(SIZE+1));
	if(gMatrixA == NULL)	return -1;
	gMatrixB = (float*)malloc(sizeof(float)*SIZE*(SIZE+1));
	if(gMatrixB == NULL)	return -1;

	// For test
	gMatrixResult = (float*)malloc(sizeof(float)*SIZE*(SIZE+1));
	if(gMatrixResult == NULL)	return -1;
	memset((void*)gMatrixResult, 0x00, sizeof(float)*SIZE*(SIZE+1));

	// For test
	gTargetRow = (float*)malloc(sizeof(float)*TILE_WIDTH);
	gCoeffi = (float*)malloc(sizeof(float)*TILE_WIDTH);

	
	srand48(time(0));
	for(i=0; i<SIZE; ++i){
		for(j=0; j<SIZE+1; ++j){
			gMatrixA[i*ROW_SIZE+j] = gMatrixB[i*ROW_SIZE+j] = drand48();
		}
	} 
	return 0;
}

int fini()
{
	free(gMatrixA);
	free(gMatrixB);

	return 0;
}

void set_option(int argc, char** argv)
{
	int opt;
	if(argc!=3 && argc!=5 ){
		printf("Usage: %s [options]\n", argv[0]);
		printf("\t-s\t--size\t\t Size of Matrix\n");
		exit(0);
	}
	do{ 
		opt = getopt_long(argc, argv,"s:t:", input_option, NULL); 
		switch(opt)
		{
			case 's':
				SIZE = atoi(optarg);
				ROW_SIZE = SIZE+1;
				break;
			case -1:
				break;
			default:
				abort();
		}
	}while(opt!=-1);
}

void print_result()
{
	printf("\nSize: %d\n", SIZE);
	printf("L2Norm: %f\n", dL2Norm);
	printf("Time Ellapsed: %f msec\n", (endTime.tv_sec - startTime.tv_sec)*1000 + (endTime.tv_usec - startTime.tv_usec)/1000.0);
}

__global__ void gaussian_elim2(float* matrixdA, float* matrixdResult, int SIZE)
{
	int ROW_SIZE = SIZE+1;
	int current;
	int tx = threadIdx.x;
	int ty = threadIdx.y;

	__shared__ float sm_target_row[TILE_WIDTH];
	__shared__ float sm_coeffi[TILE_WIDTH];

	int current_target = 0;

	// 1) Store target row value into shared memory
	if(ty == 0){
		sm_target_row[tx] = matrixdA[current_target*ROW_SIZE+tx];
	}
	__syncthreads();

	// 2) Store coeffi value into shared memory
	if(tx == 0){
		sm_coeffi[ty] = matrixdA[ty*ROW_SIZE+tx] / sm_target_row[tx];
	}
	__syncthreads();

	// 3) Compute all rows
	if(ty != current_target)
		matrixdA[ty*ROW_SIZE+tx] = matrixdA[ty*ROW_SIZE+tx] - sm_coeffi[ty]*sm_target_row[tx];

	__syncthreads();

	return;
}

__global__ void gaussian_elim3(float* matrixdA, float* matrixdResult, int SIZE, float* dCoeffi, float* dTargetRow)
{
	int ROW_SIZE = SIZE+1;
	int tx = threadIdx.x;
	int ty = threadIdx.y;

	__shared__ float sm_target_row[TILE_WIDTH];
	__shared__ float sm_coeffi[TILE_WIDTH];

	int current_target = 0;

	// 1) Store target row value into shared memory
	if(ty == 0){
		sm_target_row[tx] = matrixdA[current_target*ROW_SIZE+tx];
		dTargetRow[tx] = sm_target_row[tx];
	}
	__syncthreads();

	// 2) Store coeffi value into shared memory
	if(tx == 0){
		sm_coeffi[ty] = matrixdA[ty*ROW_SIZE+tx] / sm_target_row[tx];
		dCoeffi[ty] = sm_coeffi[ty];
	}
	__syncthreads();

	// 3) Compute all rows
	if(ty != current_target)
		matrixdA[ty*ROW_SIZE+tx] = matrixdA[ty*ROW_SIZE+tx] - sm_coeffi[ty]*sm_target_row[tx];

	__syncthreads();

	return;
}

__global__ void gaussian_elim4(float* matrixdA, float* matrixdResult, int SIZE)
{
	int current_target;
	int ROW_SIZE = SIZE+1;
	int idxCol = threadIdx.x;
	int idxRow = threadIdx.y;

	__shared__ float sm_target_row[TILE_WIDTH];
	__shared__ float sm_coeffi[TILE_WIDTH];

	for(current_target= 0; current_target<TILE_WIDTH; ++current_target){

		// 1) Store target row value into shared memory
		if(idxRow == current_target){
			sm_target_row[idxCol] = matrixdA[current_target*ROW_SIZE+idxCol];
		}
		__syncthreads();

		// 2) Store coeffi value into shared memory
		if(idxCol == current_target){
			sm_coeffi[idxRow] = matrixdA[idxRow*ROW_SIZE+idxCol] / sm_target_row[idxCol];
		}
		__syncthreads();

		// 3) Compute all rows
		if(idxRow != current_target)
			matrixdA[idxRow*ROW_SIZE+idxCol] = matrixdA[idxRow*ROW_SIZE+idxCol] - sm_coeffi[idxRow]*sm_target_row[idxCol];

		__syncthreads();
	}

	return;
}


int device_malloc()
{
	// Copy to device
	if(cudaSuccess != cudaMalloc(&dMatrixA, sizeof(float)*SIZE*(SIZE+1))){
		ERROR("Fail to cudaMalloc()");
		return -1;
	}
	if(cudaSuccess != cudaMalloc(&dMatrixResult, sizeof(float)*SIZE*(SIZE+1))){
		ERROR("Fail to cudaMalloc()");
		return -1;
	}
	if(cudaSuccess != cudaMalloc(&dTargetRow, sizeof(float)*TILE_WIDTH)){
		ERROR("Fail to cudaMalloc()");
		return -1;
	}
	if(cudaSuccess != cudaMalloc(&dCoeffi, sizeof(float)*TILE_WIDTH)){
		ERROR("Fail to cudaMalloc()");
		return -1;
	}
	cudaMemcpy(dMatrixA, gMatrixA, sizeof(float)*SIZE*(SIZE+1), cudaMemcpyHostToDevice);

	return 0;
}

int device_free()
{
	// Free device memory
	cudaFree(dMatrixA);
	cudaFree(dMatrixResult);
	cudaFree(dCoeffi);
	cudaFree(dTargetRow);
	return 0;
}


int main(int argc, char** argv)
{
	set_option(argc, argv);

	// Host Initialize
	if(init()<0){
		ERROR("Fail to init()"); 
		return -1;	
	}

	show_matrix(gMatrixA);
	gettimeofday(&startTime, NULL);

#if 1
	if(device_malloc() <0){
		ERROR("Fail to device_malloc()");
		return -1;
	}

	dim3 dimBlock(SIZE, SIZE);
	dim3 dimGrid(1, 1);

	// gaussian_elim3<<<dimGrid, dimBlock>>>(dMatrixA, dMatrixResult, SIZE, dCoeffi, dTargetRow);
	gaussian_elim4<<<dimGrid, dimBlock>>>(dMatrixA, dMatrixResult, SIZE);
	gaussian_elim4<<<dimGrid, dimBlock>>>(dMatrixA, dMatrixResult, SIZE);

	// Copy to Host
	cudaMemcpy(gMatrixA, dMatrixA, sizeof(float)*SIZE*(SIZE+1), cudaMemcpyDeviceToHost);
	cudaMemcpy(gCoeffi, dCoeffi, sizeof(float)*TILE_WIDTH, cudaMemcpyDeviceToHost);
	cudaMemcpy(gTargetRow, dTargetRow, sizeof(float)*TILE_WIDTH, cudaMemcpyDeviceToHost);

	show_matrix(gMatrixA);

	device_free();
#endif	


	// show_matrix(gMatrixA);

	gettimeofday(&endTime, NULL);

	print_result();

	// Host Finalize
	fini();
	return 0;
}
